Appendix 3: Additional models results

#Loading data
load("scripts/datasets.Rdata")
load("models_outputs/models_combined.Rdata") # saved models results combined traits
load("models_outputs/models_traits.Rdata") # saved models results separate traits
#Loading auxiliary function for R2 calculations
source("scripts/r2_auxiliary_function.R")
# variance inflation factor function
vif.mer <- function (fit) {
    ## adapted from rms::vif
    v <- vcov(fit)
    nam <- names(fixef(fit))
    ## exclude intercepts
    ns <- sum(1 * (nam == "Intercept" | nam == "(Intercept)"))
    if (ns > 0) {
        v <- v[-(1:ns), -(1:ns), drop = FALSE]
        nam <- nam[-(1:ns)]
    }
    d <- diag(v)^0.5
    v <- diag(solve(v/(d %o% d)))
    names(v) <- nam
    v
}

1. Models with traits

Specification of the models: We used lme4 package to perform a GLMM with binomial (proportion) distribution. An example of the code for each dataset are as follows:

mhigh.spe <- glmer(cbind(occor, n.visit-occor) ~  
                    forest_site400*lbody_size + 
                    forest_site400*nest +
                    forest_site400*diet + 
                    forest_site400*lower_stratum +
                    forest_land*lbody_size + 
                    forest_land*nest + 
                    forest_land*diet +
                    forest_land*lower_stratum + 
                    (forest_site400 + forest_land|sp) + 
                    (1|landscape:sp) + (1|site:sp) + 
                    (lbody_size + nest + diet + lower_stratum|landscape) +
                    (lbody_size + nest + diet + lower_stratum|site),
                    family=binomial, data=high.spe,
                    nAGQ = 1, control = glmerControl(optimizer = "bobyqa",
                                    optCtrl = list(maxfun = 500000)))

We ran separate models for each assemblage and trait. Afterwards, we ran one model with the combination of the traits body mass, diet, nest type and % of lower strata use. Table S3. shows the marginal R2 of all models terms.

traits = c("lbody_size", "nest", "diet", "lower_stratum")
rhigh.spe <- rbind(r2.calc(mhigh.spe, atrib = traits) , 
                   r2.calc(bodyhigh.spe, atrib="lbody_size"), 
                   r2.calc(nesthigh.spe, atrib="nest"),
                   r2.calc(diethigh.spe, atrib="diet"),
                   r2.calc(frugihigh.spe, atrib="frugivory"),
                   r2.calc(insehigh.spe, atrib="insectivory"),
                   r2.calc(low_strathigh.spe, atrib = "lower_stratum"),
                   r2.calc(stratumhigh.spe, atrib="stratum"))
rhigh.spe <- cbind(model = rep(c("Combined", "body mass", "nest type", "main diet", 
                          "% frugivory", "% insetivory","% lower strata", 
                          "foraging stratum")),
                   rhigh.spe)  %>%
  mutate_if(is.numeric, function(x) 100*x)
colnames(rhigh.spe) <- c("Model", "Total", "trait*env",  "env|sp", "lands:sp",
                         "site:sp", "trait|lands", "trait|site")

rlow.spe <- rbind(r2.calc(mlow.spe, atrib = traits) , 
                   r2.calc(bodylow.spe, atrib="lbody_size"), 
                   r2.calc(nestlow.spe, atrib="nest"),
                   r2.calc(dietlow.spe, atrib="diet"),
                   r2.calc(frugilow.spe, atrib="frugivory"),
                   r2.calc(inselow.spe, atrib="insectivory"),
                   r2.calc(low_stratlow.spe, atrib = "lower_stratum"),
                   r2.calc(stratumlow.spe, atrib="stratum"))
rlow.spe <- cbind(model = rep(c("Combined", "body mass", "nest type", "main diet", 
                                 "% frugivory", "% insetivory","% lower strata", 
                                 "foraging stratum")),
                   rlow.spe)  %>%
  mutate_if(is.numeric, function(x) 100*x)
colnames(rlow.spe) <- c("Model", "Total", "trait*env",  "env|sp", "lands:sp",
                         "site:sp", "trait|lands", "trait|site")
rhigh.gen <- rbind(r2.calc(mhigh.gen, atrib = traits) , 
                   r2.calc(bodyhigh.gen, atrib="lbody_size"), 
                   r2.calc(nesthigh.gen, atrib="nest"),
                   r2.calc(diethigh.gen, atrib="diet"),
                   r2.calc(frugihigh.gen, atrib="frugivory"),
                   r2.calc(insehigh.gen, atrib="insectivory"),
                   r2.calc(low_strathigh.gen, atrib = "lower_stratum"),
                   r2.calc(stratumhigh.gen, atrib="stratum"))
rhigh.gen <- cbind(model = rep(c("Combined", "body mass", "nest type", "main diet", 
                                 "% frugivory", "% insetivory","% lower strata", 
                                 "foraging stratum")),
                   rhigh.gen)  %>%
  mutate_if(is.numeric, function(x) 100*x)
colnames(rhigh.gen) <- c("Model", "Total", "trait*env",  "env|sp", "lands:sp",
                         "site:sp", "trait|lands", "trait|site")

rlow.gen <- rbind(r2.calc(mlow.gen, atrib = traits) , 
                  r2.calc(bodylow.gen, atrib="lbody_size"), 
                  r2.calc(nestlow.gen, atrib="nest"),
                  r2.calc(dietlow.gen, atrib="diet"),
                  r2.calc(frugilow.gen, atrib="frugivory"),
                  r2.calc(inselow.gen, atrib="insectivory"),
                  r2.calc(low_stratlow.gen, atrib = "lower_stratum"),
                  r2.calc(stratumlow.gen, atrib="stratum"))
rlow.gen <- cbind(model = rep(c("Combined", "body mass", "nest type", "main diet", 
                                "% frugivory", "% insetivory","% lower strata", 
                                "foraging stratum")),
                  rlow.gen)  %>%
  mutate_if(is.numeric, function(x) 100*x)
colnames(rlow.gen) <- c("Model", "Total", "trait*env",  "env|sp", "lands:sp",
                        "site:sp", "trait|lands", "trait|site")
todos <- bind_rows(list(rhigh.spe, rlow.spe, rhigh.gen, rlow.gen))

kable(todos, "latex", booktabs=T, caption= "Overall and marginal R-squared of trait models in each dataset. For the marginal R-squared terms see Table 2 (main text).") %>%
   group_rows("Specialists", 1, 16) %>%
   group_rows("     High quality", 1,8) %>%
   group_rows("     Low quality", 9,16) %>%
   group_rows("Generalists", 17, 32) %>%
   group_rows("     High quality", 17,24) %>%
   group_rows("     Low quality", 25,32)

2. Models coeficients

Tables S3., S3., S3., and S3. show the coefficients for each model.

 tidy(mhigh.spe, effects = "fixed") %>% 
  kable(digits=2,"latex", booktabs=T, caption = "Fixed effects oefficients for the model of specialists in high-quality matrix landscapes.")
 tidy(mlow.spe, effects = "fixed") %>% 
  kable(digits=2,"latex", booktabs=T,  caption = "Fixed effects oefficients for the model of specialists in low-quality matrix landscapes.")
 tidy(mhigh.gen, effects = "fixed") %>% 
  kable(digits=2, "latex", booktabs=T, caption = "Fixed effects oefficients for the model of generalists in high-quality matrix landscapes.")
 tidy(mlow.gen, effects = "fixed") %>% 
  kable(digits=2,"latex", booktabs=T, caption = "Fixed effects oefficients for the model of generalists in low-quality matrix landscapes.")

3. Models diagnostic

Variance Inflation Factor of the model parameters for each dataset in Table S3..

vifhigh.spe <- data.frame(termo = names(vif.mer(mhigh.spe2)),
                        high.spe =vif.mer(mhigh.spe2))
viflow.spe <- data.frame(termo = names(vif.mer(mlow.spe2)),
                        low.spe =vif.mer(mlow.spe2))
viflow.gen <- data.frame(termo = names(vif.mer(mhigh.gen2)),
                        low.gen =vif.mer(mhigh.gen2))
viflow.gen <- data.frame(termo = names(vif.mer(mlow.gen2)),
                        low.gen =vif.mer(mlow.gen2))
varif <- vifhigh.spe %>% full_join(viflow.spe, "termo") %>% 
  full_join(viflow.gen, "termo") %>% full_join(viflow.gen, "termo")
varif$termo <- c("forest.local", "body_mass", "nest_closed", "nest_open_semi",
                 "diet_insectivorous", "diet_onivorous", "lower_strata", 
                 "diet_granivorous", "forest.landscape", "diet_nectarivorous")
colnames(varif) <-c("parameter","Coffee", "lowture", "Coffee", "lowture")

kable(varif, "latex", booktabs=T, caption= "Variance Inflation Factor index for combined traits models in each dataset.", digits=2) %>%
 add_header_above(c(" " = 1, "Specialists" = 2, "Generalists"= 2))

Example of the residual diagnostic of the model with the combined traits (main diet, body mass, nest type and % of lower strata use) for the forest specialists in high-quality matrix landscapes. The models’ diagnostics for the other assemblages were all similar and can be checked in this Rmd file.

Residual correlations among species and sites

Below we present the Kendall correlations for the residuals among species and sites for the models using the predictions for site:sp random effect (Observation Level Random Effect). For the residual correlations we followed the code provided by Miller, Damschen & Ives (2018).

nsites = 40
nspp = length(unique(high.spe$sp))
dat <- high.spe

dat$resid <- as.matrix(ranef(mhigh.spe)$`site:sp`)
X <- matrix(dat$resid, nrow=nsites, ncol=nspp, byrow=F)
rownames(X) <- unique(dat$ponto)
colnames(X) <- unique(dat$sp)

# species correlations
corrs.sp <- cor(X, method="kendall")
for(i in 1:nspp) corrs.sp[i,i] <- NA

# site correlations
corrs.site <- cor(t(X), method="kendall")
for(i in 1:nsites) corrs.site[i,i] <- NA

Range of species correlations: -0.4, 0.43. Range of sites correlations: -0.3, 0.27.

corrplot(corrs.sp, method="square", type="upper",diag=F,tl.cex=0.5,cl.cex=0.5)
Species residual Kendall correlations for the specialist species in high-quality matrix landscapes.

Species residual Kendall correlations for the specialist species in high-quality matrix landscapes. 

corrplot(corrs.site, method="square", type="upper",diag=F,tl.cex=0.5,cl.cex=0.5)
Sites residual Kendall correlations for the specialist species in high-quality matrixlandscapes.

Sites residual Kendall correlations for the specialist species in high-quality matrixlandscapes. 

x.sp <- matrix(corrs.sp, ncol=1)
x.site <- matrix(corrs.site, ncol=1)
x.sp <- x.sp[!is.na(x.sp)]
x.site <- x.site[!is.na(x.site)]

# Figure histogram
par(mfrow=c(1,2), mai=c(1,1,.5,.2))
hist(corrs.sp, main="", probability=T, breaks=.1*(-10:10), xlab="Correlation coefficients", ylab="Probability", ylim=c(0,4), cex.lab=1.2)
mtext("Species", side=3,cex=1.2)

hist(corrs.site, main="", probability=T, breaks=.1*(-10:10), xlab="Correlation coefficients", ylab="Probability", ylim=c(0,4), cex.lab=1.2)
mtext("Sites", side=3, cex=1.2)
Histograms of the residual Kendall correlations for the specialists species in high-quality matrix landscapes.

Histograms of the residual Kendall correlations for the specialists species in high-quality matrix landscapes.

Residual diagnostic

We used DHARMa package (Hartig (2018)) for the diagnostic of quantile residuals.

resce <- simulateResiduals(mhigh.spe, n=1000, ref.form=~0)
plot(resce)

Residuals against predictors:

par(mfrow=c(4,2))
plotResiduals(resce, high.spe$forest_site400, xlab="Local forest cover 400m")
plotResiduals(resce, high.spe$forest_land, xlab="Landscape forest")
plotResiduals(resce, high.spe$lbody_size, xlab="Body size")
plotResiduals(resce, high.spe$nest, xlab="Nest type")
plotResiduals(resce, high.spe$diet, xlab="Main diet")
plotResiduals(resce, high.spe$low_strat, xlab="% lower strata use")

Predictions for each species local forest cover

Landscape forest cover was fixed in 30%.

#landscape cover 30%
fland = (30 - mean(high.spe$forest_landOrig))/
    sd(high.spe$forest_landOrig)
# local forest cover
f400st.high <- seq(min(high.spe$forest_site400), max(high.spe$forest_site400),
                     length.out=5)
f400Orig.high <- f400st.high * sd(high.spe$forest_site400Orig) +
                                              mean(high.spe$forest_site400Orig)
local.high <- data.frame(forest_site400=f400st.high, 
                         forest_site400Orig = f400Orig.high)

f400st.low <- seq(min(low.spe$forest_site400), max(low.spe$forest_site400),
                    length.out=5)
f400Orig.low <-f400st.low * sd(low.spe$forest_site400Orig) +
                  mean(low.spe$forest_site400Orig)
local.low <- data.frame(forest_site400=f400st.low, 
                         forest_site400Orig = f400Orig.low)
# high quality specialists
sp.tra <- high.spe %>% select(sp, lbody_size, nest, diet, lower_stratum) %>%
  distinct() 
shigh.spe <- expand.grid(sp = unique(high.spe$sp),
                    forest_land = fland,
                    forest_site400 =f400st.high,
                    landscape = "P02",
                    site = "P02.P00") %>%
  left_join(sp.tra, "sp") %>%
  left_join(local.high, "forest_site400")
shigh.spe$pretes <- predict(mhigh.spe, newdata=shigh.spe, type="response",
                       re.form=NULL)

# low quality specialists
sp.tra2 <- low.spe %>% select(sp, lbody_size, nest, diet, lower_stratum) %>%
  distinct()
slow.spe <- expand.grid(sp = unique(low.spe$sp),
                    forest_land = fland,
                    forest_site400 =f400st.low,
                    #forest_site400Orig = f400Orig.low,
                    landscape = "148",
                    site = "148.P10") %>%
  left_join(sp.tra2, "sp")%>%
  left_join(local.low, "forest_site400")
slow.spe$pretes <- predict(mlow.spe, newdata=slow.spe, type="response",
                           re.form=NULL)
bind_rows(list(High=shigh.spe,Low=slow.spe), .id="Matrix quality") %>%
  mutate(forest_site400Orig = round(forest_site400Orig)) %>%
  ggplot(aes(x=forest_site400Orig, y=pretes, col=`Matrix quality` )) +
  geom_line() +
  facet_wrap(~sp, ncol=8) + 
  theme_cowplot() +
  theme(legend.position = "top",
        strip.text.x = element_text(hjust=0, size=10))+
 scale_y_continuous(name= "Occurrence probability") +
  scale_x_continuous(name="Local forest cover (%)") +
  ggtitle("Forest specialists")

# high quality generalists
sp.tra <- high.gen %>% select(sp, lbody_size, nest, diet, lower_stratum) %>%
  distinct() 
shigh.gen <- expand.grid(sp = unique(high.gen$sp),
                    forest_land = fland,
                    forest_site400 =f400st.high,
                    landscape = "P02",
                    site = "P02.P00") %>%
  left_join(sp.tra, "sp") %>%
  left_join(local.high, "forest_site400")
shigh.gen$pretes <- predict(mhigh.gen, newdata=shigh.gen, type="response",
                       re.form=NULL)

# low quality  generalists
sp.tra2 <- low.gen %>% select(sp, lbody_size, nest, diet, lower_stratum) %>%
  distinct()
slow.gen <- expand.grid(sp = unique(low.gen$sp),
                    forest_land = fland,
                    forest_site400 =f400st.low,
                    #forest_site400Orig = f400Orig.low,
                    landscape = "148",
                    site = "148.P10") %>%
  left_join(sp.tra2, "sp")%>%
  left_join(local.low, "forest_site400")
slow.gen$pretes <- predict(mlow.gen, newdata=slow.gen, type="response",
                           re.form=NULL)
bind_rows(list(High=shigh.gen,Low=slow.gen), .id="Matrix quality") %>%
  mutate(forest_site400Orig = round(forest_site400Orig)) %>%
  ggplot(aes(x=forest_site400Orig, y=pretes, col=`Matrix quality` )) +
  geom_line() +
  facet_wrap(~sp, ncol=8) + 
  theme_cowplot() +
  theme(legend.position = "top",
        strip.text.x = element_text(hjust=0, size=10))+
 scale_y_continuous(name= "Occurrence probability") +
  scale_x_continuous(name="Local forest cover (%)") +
  ggtitle("Forest generalists")

References

Hartig, F. (2018). DHARMa: Residual Diagnostics for Hierarchical (Multi-Level / Mixed) Regression Models.
Miller, J.E.D., Damschen, E.I. & Ives, A.R. (2018). Functional traits and community composition: A comparison among community-weighted means, weighted correlations, and multilevel models. Methods in Ecology and Evolution 0.